#!/usr/bin/perl -w
use strict;
use Term::ANSIColor;


#+=============================================================================+
# FILE:     matchRNA.pl
# SUMMARY:  Illustrate where results of RNAfold or CentroidFold lie 
#			compared to full sequence
# AUTHOR:   Ravi Mehra
# CREATED:  Jan. 20, 2010
# MODIFIED:
# NOTE:     Ver.1
#
# Usage: matchRNA.pl file1.pri file2.01v100 file3.01v100m8 
#
#   INPUT 1: file1.pri (fasta formated)
#   INPUT 2: file2.01v100 (blast matches)
#
#   OUTPUT 
# 
#+=============================================================================+
#        1         2         3         4         5         6         7         8
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#+=============================================================================+


#-------------------------------------------------------------------------------
use Cwd;
use Getopt::Long;
#-------------------------------------------------------------------------------

my ($format, $help, $minHitCount);

Getopt::Long::Configure ('bundling');
GetOptions ('f' => \$format, "help|?|h" => \$help, 'c|count=i' => \$minHitCount); 

if (!defined($minHitCount)) {
	$minHitCount = 0;
}

# Check to make sure all three files are given
if (@ARGV < 2 or $help) {
  print " Usage: \"./matchRNA.pl file1.pri file2.01v100\"\n";
  print " Options -f : Format for file output, i.e. \"./matchRNA.pl -f pri_file blast_file > matches.txt\"\n";
  print " Maximizing the command line window is recommended.\n";
  exit(0);
}

# Open all three files for reading
my $pri_file = $ARGV[0];
open (PRI,  "<$pri_file")  || die "Cannot open $pri_file!\n";
my $blast_file = $ARGV[1];
open (BLAST,  "<$blast_file")  || die "Cannot open $blast_file!\n";
#my $table_file = $ARGV[2];
#open (TABLE,  "<$table_file")  || die "Cannot open $table_file!\n";

# The various hashes used
my %subjectSeqHash = ();			# Stores the Query ID and it's associated sequence from the .pri file
my %matchHash = ();					# Stores the Subject with all of it's necessary information into an array
my %queryIDSubjectHash = ();		# Stores the Query ID and it's associated Subject
my %queriesCombinedHash = ();		# Stores the Query ID with all of it's Subject matches
my %centroidOutputHash = ();		# Stores the Query ID and it's assoicated fold shape reported by centroid fold

# Parse through the pri file extracting info for the %subjectSeqHash
print "Parsing through the .pri file ... " unless $format;
while (<PRI>) {
    # Read line 
    my ($line) = $_;
    chomp ($line);
    # Check the line begins with fasta tag format such as ">male_1_1"
    if ($line =~ /^>/) {
        my $queryID = $line;
        $queryID = substr ($queryID, 1, length($queryID));
        # Read the next line (the Sequence ID's associated sequence) if eof hasn't been reached
        my $fullSequence = $_ if defined($_=<PRI>); 
        chomp ($fullSequence);
        # Store the pair into the hash
        $subjectSeqHash {$queryID} = $fullSequence;
    }
    else {
		# This should not be encountered as there should be no other lines
        print "\nUnexpected line in .pri file: ", $line;
    }
}
print "done\n" unless $format;

# Parse through the blast matches file extracting info for the %matchHash
print "Parsing through the blast output file ... " unless $format;
my $curQueryID;
while (<BLAST>) {
    # Read line 
    my ($line) = $_;
    chomp ($line);
    
    if ( $line =~ /^Query=/ ) {
		my @queryIDLineArray = split(' ', $line);
		$curQueryID = $queryIDLineArray[1];
	}
    # Check the line begins with fasta tag format such as ">male_1_1"
    if ( $line =~ /^>/ ) {
        my $subject = $line;
        my $subjectPercent = 0;
        # Get the subject ID
        $subject = substr ($subject, 1, length($subject));
        my $subjectCount = substr($subject, rindex($subject, "_") + 1, length($subject));	# Gets the count of the subject
        my $subjectShort = substr($subject, 0, index($subject, "_")); # Gets just the first part of the name
        # Since there are other lines in the file we are not interested in, it will skip them
        while (<BLAST>) {
            my ($subLine) = $_;
            chomp ($subLine);
            if ($subLine =~ /Identities =/) {
				# Get just the percentage from the line
				$subjectPercent = substr($subLine, rindex($subLine, "(") + 1);
				$subjectPercent = substr($subjectPercent, 0, length($subjectPercent)-1);
			}
            # Skips everything until it reaches the "Query:" line
            if ($subLine =~ /^Query:/) {
				# Extract the query line's info
                my @queryLineArray = split(' ', $subLine);
                my $querySeq = $queryLineArray[2];
                my $querySeqStart = scalar($queryLineArray[1]);
                my $querySeqEnd = scalar($queryLineArray[3]);
                # Get the next line if eof is not reached
                my $matches = $_ if defined($_=<BLAST>);
                chomp ($matches);
                # Stores the hash marks which indicate nt matches
                $matches =~ s/\s+//;
                # Get the next line if eof is not reached
                my $subjectSeqLine = $_ if defined($_=<BLAST>);
                my @subjectSeqLineArray = split(' ', $subjectSeqLine);
                my $subjectSeq = $subjectSeqLineArray[2];
                my $subjectSeqStart = scalar($subjectSeqLineArray[1]);
                my $subjectSeqEnd = scalar($subjectSeqLineArray[3]);
                $matchHash {$subject} = [$querySeq, $querySeqStart, $querySeqEnd,
											$matches, uc($subjectSeq), $subjectSeqStart, 
											$subjectSeqEnd, $subjectCount, $subjectShort, $subjectPercent];
				$queryIDSubjectHash {$subject} = $curQueryID;
                last;
            }
        }
    }
    else {
        # Other miscelaneous lines are ignored
    }
}
print "done\n" unless $format;

=pod
# read input file and main procedure
print "Parsing through the m8 table formated file ... ";
while (<TABLE>) {
    # read line 
    my ($line) = $_;
    chomp ($line);
	my @lineArray = split(' ', $line);
	my $queryID = $lineArray[0];
	my $subject = $lineArray[1];
	#print $queryID, " : ", $subject, "\n";
	$queryIDSubjectHash {$subject} = $queryID;
}
print "done\n";
=cut

for my $subject (sort keys %matchHash) {
	my $queryID = $queryIDSubjectHash{$subject};
	next if (${ $matchHash{$subject} }[7] < $minHitCount);
	#print "$subject => $queryID, ";
	my @matchArray = @{ $matchHash{$subject} };
	#print ("x", $matchArray[7], "\n\n");
	my @subjectIDs = @{ $queriesCombinedHash {$queryID} } if ($#{ $queriesCombinedHash {$queryID} });
	my $i = 0;
	my $inserted = 0;
	foreach my $subjectID (@subjectIDs) {
		if ($matchArray[7] < ${ $matchHash{$subjectID} }[7]) {
			splice @{ $queriesCombinedHash {$queryID} }, $i, 0, $subject;
			$inserted = 1;
			last;
		}
		$i++ if (${ $matchHash{$subjectID} }[7] < $minHitCount);
	}
	if ($inserted == 0) {
		push @{ $queriesCombinedHash {$queryID} }, $subject; 
	}
}


print "Getting RNA structures from CentroidFold ... " unless $format;
my $centroidOutput = `centroid_fold $pri_file`;
my @centroidOutputArray = split("\n", $centroidOutput);
foreach (my $i = 0; $i < ($#centroidOutputArray + 1); $i += 3) {
    # Check the line begins with fasta tag format such as ">male_1_1"
    if ($centroidOutputArray[$i] =~ /^>/) {
        my $queryID = $centroidOutputArray[$i];
        $queryID = substr ($queryID, 1, length($queryID));
        my $foldingOutput = $centroidOutputArray[$i + 2];
        $foldingOutput = substr ($foldingOutput, 0, length($centroidOutputArray[$i + 1]));
        $centroidOutputHash {$queryID} = $foldingOutput;
    }
    else {
		# This should not be encountered as there should be no other lines
        print "\nUnexpected line in centroid output: ", $centroidOutputArray[$i];
    }
}
print "done\n\n" unless $format;


for my $queryID (sort keys %queriesCombinedHash) {
	my @subjectIDs = @{ $queriesCombinedHash{$queryID} };
	#print "$queryID => @subjectIDs\n"; 
	print ">", $queryID, "\n";
	foreach my $subjectID (@subjectIDs) {
		#print $subjectID, "\n";
		my ($querySeq, $querySeqStart, $querySeqEnd, $matches,
				$subjectSeq, $subjectSeqStart, $subjectSeqEnd, 
				$subjectCount, $subjectShort, $subjectPercent)  = @{ $matchHash{$subjectID} };
		foreach (my $i = 1; $i < ($querySeqStart - length($subjectSeqStart)); $i++) {
			print " ";
		}
		if ($subjectPercent ne "100%") {
			print $subjectSeqStart;
			my $offset = 0;
			my $lastResult = 0;
			my $result = index($matches, ' ', $offset);
			while ($result != -1) {
				print substr($subjectSeq, $lastResult, $result);
				print color 'bold red' unless $format;
				print substr($subjectSeq, $result, 1);
				print color 'reset' unless $format;
				$offset = $result + 1;
				$lastResult = $result;
				$result = index($matches, ' ', $offset);
			}
			print substr($subjectSeq, $lastResult + 1, length($subjectSeq));
			print "$subjectSeqEnd, $subjectShort, x$subjectCount, $subjectPercent\n";
		}
		else {
			print "$subjectSeqStart$subjectSeq$subjectSeqEnd, $subjectShort, x$subjectCount, $subjectPercent\n";
		}
	}
	print $subjectSeqHash{$queryID}, "\n";
	print $centroidOutputHash{$queryID}, "\n\n";
}
print "\n";

# close open files
close(PRI);
close(BLAST);
#close(TABLE);

# done
